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Introducing the effect of extinction into tlie so-called repli- 
cator equations in mathematical biology, we construct a gen- 
eral model of ecosystems. The present model shows mass 
extinction by its own extinction dynamics when the system 
initially has a large number of species (diversity). The ex- 
tinction dynamics shows several significant features such as 
a power law in basin size distribution, induction time, etc. 
The present theory can be a mathematical foundation of the 
species-area effect in the paleontologic theory for mass extinc- 
tion. 

PACS numbers: 87.10. -fe, 05.40. -Hj, 64.60.Lx 



Mechanisms of mass extinction of species in ecosys- 
tems have been studied by a number of researchers [|l],p| . 
Their conclusions can be divided into two categories, one 
emphasizing exogenous shocks |^-^ and the other, en- 
dogenous causes |3^- Building on both views, we con- 
struct a general mathematical biological model of ecosys- 
tems. This model reflects the former view, e.g., the sit- 
uation where several biotas which have been separated 
from each other for a long time are suddenly integrated 
into a larger ecological network (biotic fusion) by some 
exogenous shock |l^,|ll|. We assume that the interac- 
tion coefficients for this newly produced ecosystem can 
be written in the form of a random matrix |l4-p^. Also, 
following the latter view, we adopt the concept of an 
'extinction threshold', which we introduce into the repli- 
cator equations |16| of population dynamics. Using this 
replicator equation model with random interaction and 
the extinction threshold, we find several significant new 
features characterizing the nature of mass extinction. 

We investigate the following Nj dimensional ordinary 
differential equations: 
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These equations, the replicator equations, are generally 
used to describe the evolution of self-replicating entities, 



so-called replicators ||l^. The equations are also termed 
the game dynamical equations in the game theory. More- 
over, they are equivalent to the Ni — 1 dimensional Lotka- 
Volterra equations [|l6|. Therefore, we will use the term 
species when referring to replicators hereafter. The vari- 
able Xi denotes the population density of the species i. Nj 
denotes the initial number of species, that is, the initial 
value of the diversity. The (i,j)-th element of the matrix 
A — [aij) determines the effect of species j on the growth 
rate of species i. Here we use an = —1 for intraspecies 
interaction coefficients and assign the interspecies ones 
Oij (i 7^ j) as time-independent Gaussian random num- 
bers with mean and variance v. This assumption of ran- 
dom interactions is based on the hypothesis that a biotic 
fusion reorganizes species relationships in a random fash- 
ion pT[ . This kind of ecosystem with random interaction 
also can be produced by the so-called species-area effect 
that paleontologists have asserted to be a trigger of mass 
extinction ||l8|. For example, because the species-area 
effect may be caused by declining sea levels, which con- 
fines many biologically-isolated species to a narrow area, 
drives them into competition and, eventually, brings bi- 
otic fusion. We also believe that the random interaction 
model is important as a first step of understanding the 
behavior of a large ecosystem with many species inter- 
acting in a complex way. For such a random interac- 
tion model, the local stability condition has been clarified 
||l2|-[l4| in the TV/ — > oo limit. However, the global behav- 
ior is hardly treated analytically, because the equations 
are highly nonlinear and the dynamics often shows com- 
plex behavior, such as heteroclinic orbits |19| , p0| or chaos 
pl| , p2| , even at a small degree of freedom (iV/ > 4) . 

Here we should note that extinction is not well-defined 
in the above model with a large Nj because of the hete- 
roclinic orbits. The reason is that even though a hetero- 
clinic orbit approaches a saddle where some species are 
extinct, the population densities never reach exactly zero 
and the orbit eventually leaves for another saddle where 
the population densities revive. Moreover, this transition 
among saddles continues cyclically or chaotically. In this 
sense, heteroclinic orbits have never been believed to be 
biologically siginificant. 

Considering the above problem, we introduce a param- 
eter (5, the extinction threshold, to the dynamics (0)-(0) 
: at each discrete time step, the population density Xk is 



set to zero if this quantity becomes less than i5. The pop- 
ulation densities of the surviving species {xi} {i ^ k) are 
then renormalized to satisfy X^i^^fc ^« ~ 1- This renor- 
malization implies that the niche of an extinct species 
is divided among the survivors. The diversity decreases 
through the above process, and we denote its value by N. 
Thus, the present model can be interpreted as a dynam- 
ical system whose degree of freedom is a time-dependent 
variable. Although this time-dependent degree of free- 
dom is inevitable not only in population dynamics [ p3|j24| 
but also in many fields, such a highly non-linear model 
has never been systematically analyzed. 
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FIG. 1. Basin-size distributions for (a) Nj = 64 sampled 
from 100000 initial states and averaged over 10 samples of A, 
(b) Ni = 128 from 20000 initial states and 3 samples of A. 

Whenever given a set of parameters A, S and the initial 
diversity Nj = N{t — 0), the initial state {a;i(0)} evolves 
until a steady state is achieved. In this state extinction no 
longer occurs, and there remains a stable subecosystem 
with a comparatively small number of surviving species 
{core species) Np (< Nj). In this state, although almost 
all orbits converge to an equilibrium point, we also find 
periodic orbits. Chaotic orbits are very rare. Heteroclinic 
orbits are never achieved because the existence of the 
finite 6 prohibits any orbit from approaching a saddle. 

This introduction of d also introduces a finite size effect 
into the replicator equations, because 6 coincides with a 
minimum unit of reproduction of each species, and its 
reciprocal 1/6 corresponds to the permissible population 
size of the ecosystem. Let us refer to this kind of dynam- 
ics as extinction dynamics (ED). By the series of exten- 
sive numerical simulations, we investigate novel features 
of ED, especially the dependence of ED on three param- 
eters: Nj, V and S. 

From the view point on the random system theory, 
it is important to observe a typical behavior of ED by 
executing random average of quantities over samples of 
a random matrix A. Hereafter, we will in general write 
this average as {...)a- 

The first result of this Letter concerns basin-size dis- 
tribution of ED which has a large number of basins of at- 
traction. Here, we identify each 'attractor' only by com- 



position of core species, not by its trajectory. In other 
words, even if there coexist more than one isolated attrac- 
tor in a system of core species, we do not discriminate 
these attractors and we regard them to be in one basin of 
'attraction'. The reason is that in ED such coexistence 
is rare and this classification of basins of attraction also 
agrees with a classification of subecosystems appearing 
by ED. 

In order to obtain the basin-size distribution, we (a) 
iterate ED starting from a sufficient number of random 
initial states in a system with same parameters and a 
same random matrix A, (b) count basin size as the num- 
ber of initial states which converge to each 'attractor', 
and (c) make a rank-size distribution S'(n), where the 
natural number n denotes the rank of each basin and 
can reach the total number of 'attractors' found in the 
simulation. Moreover, the above process is iterated for a 
sufficient number of random matrices A with same v and 
we finally obtain a basin-size distribution {S{n))A for a 
parameter set. (S'(r7,)}^'s for various parameter sets are 
shown in Fig. 1. 

We can see the significant characteristic that {S{n))A 
follows the power law. This indicates that the phase 
space of ED is divided in such a way that the size of 
each basin of attraction resembles each term of a geo- 
metric series. Moreover, each exponent of the power de- 
pends only on Nj, not on 6 nor v. The independence 
of 6 strongly suggests that the basin-size distribution of 
the original replicator equations (ED in the limit 6 —>■ 0) 
also follows the power law. This conjecture is relevant to 
the hierarchal coexistence of infinitely many attractors in 
the replicator equations [E5[ . The power law of rank-size 
relationship with exponent near unity is often referred to 
as Zipf's law p3|, known in linguistics and diverse fields. 
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FIG. 2. Extinction curves for various values of Ni with 
5 = 0.0001 and v = 2.0. Each curve represents an average 
taken over 1000 samples of A. 

Figure 2 shows the second result of this Letter: aver- 
age diversity as a function of time, {N{t))A (the extinc- 
tion curve). Two significant characteristics can be seen 
from this figure. First, the average final diversity {Np)a 
is independent of Nj. This result implies that no matter 



how large the diversity of initial species, the average di- 
versity of species in the final state is small in comparison. 
That is Np -C iV/. In other words, when a large random 
ecosystem emerges as a result of biotic fusion, a mass 
extinction of 'size' Nj — {Np)a will occurs. Secondly, 
the avalanche of mass extinction begins after some in- 
duction time p7| tj, and ends in each case at nearly the 
same time t/j ~ 10'^(> i/). As Nj becomes larger, tj 
also becomes larger and approaches tn. Therefore, for 
sufficiently large Nj, the extinction curve shows a sharp 
drop at tj. Such an abrupt mass extinction occurring 
on a short time scale is highly relevant to the notion of 
punctuated equilibria Eq]. 




FIG. 3. (a) Extinction curves for various values of v with 
Ni = 64 and S = 0.0001, averaged over 1000 samples of A. 
(b) The average diversity of core species {Nf)a as a function 
of V. In the limit v —> 0, {Nf)a — > Ni, because this is the 
limit of no interspecies interaction(aij — 0) , and in this limit 
the right hand side of Equation (hi) becomes 0. 




FIG. 4. Extinction curves for various values of 5 with 
Ni = 64 and v = 2.0, averaged over 1000 samples of A. 

Figure 3 concerns variation of extinction curves with 
V. As V becomes larger, the induction time tj becomes 
shorter (Fig. 3a), and {Nf)a becomes smaller (Fig. 3b). 
Consequently, when the order of the interspecies inter- 



action coefficients becomes large compared with the ab- 
solute value of the intraspecies ones {{an = —1}), the 
avalanche of mass extinction begins earlier, and a smaller 
diversity of species survives. Extinction curves for sev- 
eral values of 5 are also shown in Fig. 4. From figures 2, 
3 and 4, we can conclude that {Np)a depends only on 
V, but not on Nj nor S, which contrasts the parameter 
dependence of {S{n))A only on TV/. 




Time 

FIG. 5. Time development of {f{t))A over 1000 samples of 
A with Ni = 64, w = 2.0 and S = 10"'^. Two samples of 
average fitness f{t) are also depicted as dotted curves 

Here we also mention the time development of the av- 



erage fitness fit) = J2i=i J2j=i 0'ijXi(t)xj{t) and its ran- 
dom average (/(i))^- They are depicted in Fig. 5. It 
should be noted that the average fitness takes on pos- 
itive values, except during the short period in the be- 
ginning. The final value of / ^ 0.4 is higher than the 
standard deviation a ^ 0.16 of average fitness which is 
expected for a randomly generated ecosystem with same 
diversity {Np ^ 8). Thus, more stable ecosystems are 
self-organized by ED. We also observe that (/(<)) a does 
not show monotonic increase and reaches maximum value 
at the time near tj. It suggests that, in general, the av- 
erage fitness shoots up by the avalanche of extinction of 
low fitness species around the induction time and set- 
tles down to the final value via competition among core 
species. 
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FIG. 6. Snap shots of interaction coefficients distribution 
averaged over 2000 samples of A with Nj = 64, v — 2.0 and 
(5 = 10~^ 

Time development of distribution of interspecies inter- 
action coefficients aij{i ^ j) among surviving species is 
depicted in Fig. 6. The average of Uij shifts to positive 
value along the time, which means that the relationship 
among the species becomes cooperative by ED. This also 
contributes the increase of average fitness. It should be 
noted that the distribution continuously holds its shape 
of gauss distribution. Therefore, the interspecies inter- 
action coefficients of core species are still random, and 
various types of relationship among core species are re- 
alized by ED. 

In this Letter, we ignore any effects by immigrants 
or invaders, which increase diversity, and we focus on 
global biotic fusion where no species ever comes from out- 
side. Moreover, we do not consider any mutants because 
avalanche of mass extinction occurs so quickly that any 
evolution of mutants never follows. By neglecting these 
effects, the nature of extinction on rather short time scale 
is exclusively clarified. However, by introducing the ef- 
fect of the increasing diversity, we can study the nature 
of ED on much larger time scale. For example, it must be 
interesting problem whether ED shows the self-organized 
criticality [||j|9|. 

The present theory suggests that a biotic fusion by 
some external shock will cause a mass extinction if the 
fusion occurs in large scale and the interspecies interac- 
tions become random. It also can be a mathematical 
biological foundation of the species-area effect because it 
possibly plays a similar role like the biotic fusion. Fur- 
thermore, our results clarify the importance of the extinc- 
tion threshold, i.e., the finite size effect for the replicator 
equations. 

Finally, we strongly believe in the universality of our 
model, consisting of replicator equations and the finite 
size effect, because the former are accepted models in 
many scientific fields, such as sociobiology, prebiotic evo- 
lution of macromolccules, mathematical ecology, popula- 
tion genetics, game theory and even economics, and the 
existence of the latter is inevitable for such fields. 
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